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We consider a van der Pol-Mathieu (vdPM) equation with parametric forcing, which arises in a 
simplified model of duty plasma with dust-charge fluctuation We make a detailed numerical 
investigation and show that the system can be driven to chaos either through a period doubling 
cascade or though a subcritical pitchfork bifurcation over an wide range of parameter space. We also 
discuss the frequency entrainment or frequency-locked phase of the dust-charge fluctuation dynamics 
and show that the system exhibits 2:1 parametric resonance away from the chaotic regime. 
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bX). The subject of parametric excitation can be traced back to Faraday in 1831 Q, when he observed that surface 
waves in a fluid-filled cylinder under vertical excitation exhibited twice the period of the excitation itself. The most 
simplified version of parametric excitation was given by Mathieu in 1868 [3j related to the vibrations of an elliptical 
membrane, which now has become a model equation for response of many systems to sinusoidal parametric excitation. 
The simplest Mathieu equation can be stated as 

x + (5 + e cost)x = 0, (1) 
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^■"^ , with S and e as constants. There have been extensive investigations of parametric excitation and resonance related 
Mathieu equation by several authors @, 0] ■ One of the very common examples of parametric forcing modeled by 
the nonlinear Mathieu equation is the forced and unforced inverted pendulum 0]. 

In this work, we have studied the parametric excitation and resonance of the van der Pol-Mathieu (vdPM) equation, 
which arises in a simplified model of dusty plasma with dust-charge fluctuation. Dusty plasmas are characterized by 
presence of massive dust (impurities) particles embedded in an electron-ion plasma [9j. Immersed in the plasma, 
^ ' the dust particles acquire charges by collecting the electrons and ions on their surfaces, which are mostly negative. 
However, the charge on a dust particle is never a constant and varies temporally. Thus, along with other their 
dynamical properties, the charge of the dust particles becomes a dynamic variable, which can severely modify the 
plasma properties. The presence of dust particles in a plasma may modify the dynamics of the plasma in many ways, 
of which the most prominent is the appearance of low frequency dust-acoustic waves [9j . The dust-charge fluctuation 
is known to damp the acoustic waves in a dusty plasma d, |T(|. Recently, Momeni, Kourakis, and Shukla [III nas 
studied a simplified model of nonlinear dust-charge fluctuation based on a vdPM equation, where they have discussed 
the stability regions of the vdPM equation and have shown the existence of stable and unstable periodic orbits in 
different parameter space. This vdPM equation is originally derived by Saitou and Honzawa [l|, where they have 
shown that in a very restricted region of parameter space, the system exhibit chaotic behavior. They argue that the 
oscillations leading to chaos basically stems out of the balancing between the van der Pol (vdP) and Mathieu-like 
; 1 ■ terms. 

We report, in this paper, a detailed investigation of the vdPM equation for dust-charge fluctuation and show that 
the vdPM equation proposed by Saitou [l| can exhibit chaotic behavior over an wide range of relevant parameter 
space, which is, in many instances, preceded by period doubling cascades. We have found that both period doubling 
and pitchfork bifurcations take place as one varies the bifurcation parameters, both leading to chaos. The range of 
parameters for a chaotic regime comes out to be not as restrictive as pointed out by Saitou [l|. We have discussed 
the stability and bifurcation of the nonlinear system and found that the system can be completely deterministic in 
between chaotic regimes. The paper is organized as follows. In Section II, we formulate the nonlinear dust-charge 
fluctuation model yielding the vdPM equation and discuss about the parametric forcing. In Section III, we have 
discussed about frequency entrainment (frequency-locked phase) of the nonlinear oscillator, driven by the parametric 
forcing term. We have shown that far away from the chaotic regime, the system can be driven by a 2:1 parametric 
resonance leading to a stable limit cycle. Away from the resonance, it displays quasi-periodic behavior. We discuss 
the stability and bifurcation of the periodic orbit of the nonlinear system in Section IV with the help of Floquet 
stability theory and show that the instability of a limit cycle may manifest through a period doubling bifurcation. 
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In Section V, we explore the chaotic regime of the vdPM equation, where we show that the system can be driven to 
chaos over an wide range of parameters. We have found that the route to chaos may be through a period doubling 
cascade. We draw the conclusions in Section VI. 



II. DUST-CHARGE FLUCTUATION MODEL 



We consider an unmagnctizcd collisionless plasma consisting of electrons, ions, and massive dust particles, which 
become charged by acquiring charged particles (ion or electrons) on their surfaces. The subject of dust-charge 
fluctuation in a dusty plasma is a well studied process which, primarily, has a damping effect on the acoustic waves 
Here, we assume that the number densities of the charged particles are considerably larger than that of the 
dust particles, so that the effect of dust-charge fluctuation on the dynamics of the electrons and ions is negligible and 
charge neutrality is always satisfied So, the charge on the dust grains q(t) becomes a time dependent function. 

Assuming the equilibrium (unperturbed) state to be static, we can write the nonlinear continuity, momentum 
balance, and the Poisson equation for dust-charge fluctuation as [Tl. [TTI] . 

dn „ 1 „ , . . 

— +n \7-v = an--l3n i 1 2 
at 6 

dv . . 

md ~dt = ^ ' 

e V-E = qn, (4) 

where mj, uq, Eq are dust mass, equilibrium dust number density, and permittivity of free space. The variables n, v, E 
are perturbed dust density, velocity, and electric field. In the first equation, Eq. ([2]) , the terms on the right hand side 
denote the rate of production and loss of charged dust grains, where a and j3 are constants of proportionality. In 
writing these terms, we have assumed that the production rate of charged dust particles is proportional to the dust 
density. The cubic loss term appears mainly due to the loss of dust grains through a three-body recombination process 
[l|. In the momentum equation, Eq.Q, we have assumed the dust particles to be cold which basically eliminates 
any variant of dust-acoustic waves. In writing these equations, we have assumed that the average dust velocity, v, is 
fairly uniform in space and its spatial gradient is considerably smaller so that the convective derivative term in the 
momentum equation, namely, the term (v ■ V)u can be neglected. Further, we have approximated the term V • (nv) 
with noV • v, assuming a uniform distribution of the charged dust particles in space (Vn ~ OWllJ. We assume the 
dust-charge q(t) to be changing harmonically with time with a frequency v and use the ansatz [ll. 

q{t) = q (l - eA cosvt) 1/2 , (5) 

where the term (eA) denotes the strength of charge fluctuation. Note that, in principle, there is no need to restrict 
the value of the term (eA) to a smaller value, which gives us freedom to explore an wider parameter space of a-e. 
Without loss of generality, we consider only one dimension, z and write Eqs.([2]|4|) as, 

dn dv z 1 3 

-+n — = an--/**, (6) 

md ~dT = ^ ' 

£o-^— = qn. (8) 

az 

By taking a z-derivative of Eq.Q, we can eliminate the terms involving perturbed velocity and electric field using 
Eqs.{8]) and ([6]) to get a coupled differential equation in perturbed density, 

cPn , 9n dn n , 

— 2 - (a -fin 2 )— +nwj(l - eA cosirt) = 0, (9) 

(IL (lib 

where uid = ('T-o9o/ m d e o) 1 ^ 2 is the plasma frequency corresponding to the dust particles. The above equation, Eq.Q 
can be classified as van der Pol-Mathieu (vdPM) equation owing to the nonlinear term (a — j3n 2 ) which is like a 
van der Pol (vdP) term [4| and parametric forcing term (1 — eA cos ft) which like the parametric term of a classical 
Mathieu equation Q. 
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A. Parametric forcing 

As is well known from the theory of classical Mathieu equation, the parametric forcing term in Eq.© makes the 
dynamics of the dust-charge fluctuation prone to chaos [1, [T3, 0, EH- As the vdP equation has a stable limit cycle, 
we can see that Eq.© should show vdP-type behavior for large a. However the parametric forcing term may still 
drive the system unstable. In all probability, we expect the onset of chaotic behavior as e increases, which should be 
more pronounced when a -C 1. 

In absence of the parametric forcing term, the stable limit cycle of the vdP equation has frequency of 1. Another 
well known result from the analysis of Mathieu equation Q is that the origin becomes unstable when the parametric 
forcing frequency v is close to twice the frequency of the unforced oscillator. So, when both the vdP and Mathieu 
terms are present, as in Eq.©, we expect a frequency entrainment at 2:1 [l5( and the system represented by Eq.Q 
must exhibit some sort of quasi-periodic and frequency-locked (entrainment) behaviors in the parameter space of e-a 
before it can be driven to chaos. The route to chaos should be through a series of quasi-periodic regime or through a 
period-doubling cascade rather than the other universal route i.e. through intermittancy [lj|. 



III. FREQUENCY ENTRAINMENT 

Entrainment dynamics plays an important part in design engineering and many other dynamical systems [l7T ]. 
Recently, frequency entrainment is shown to exist in nonautonoums chaotic oscillators fl8j |. In this work, we consider 
the possible entrainment by the parametric term, which can lead to quasi-periodicity and finally to chaos. We consider 
the following dynamical equation for this dust-charge dynamics, 

x — (a — (3x 2 )x + — eA cosvt) = 0, (10) 

where we have replaced the dust density n by the variable x. In order to facilitate the multiple scales in the problem 
and entrainment, we assume that a ~ (3 = i « 1, a small number. As the entrainment is possible only when the 
system as far away from the chaotic regime, when e is small, we assume that e < 1. We further re-scale the time by 
t — ► ujdt and write Eq. (TT0|) as 

'x — fie(l — x 2 )x + x(l — eA cos2u;t) = 0, (11) 

where [it = 5/uJd- The strength of the parametric forcing term is given by eA with e < 1 and v — 2uj. We expect that 
the parametric forcing should result in a 2:1 subharmonic resonance, when the parametric frequency v is close to 2 or 
u> ~ 1. Note that in absence of the parametric forcing term (eA = 0), the natural frequency of the oscillator is unity. 

We now introduce two different time-scales [HI, [l6| , the stretched time £ = tot and the slow time r/ — et and expand 
the forcing frequency ui about the natural frequency of the oscillator i.e. 1 with e as the expansion parameter, 

lu = l + ke + 0{e 2 ), (12) 

where A: is a detuning parameter at order e. The variable x now is expanded in a power series 

x = x (£,r 1 )+ex 1 (Z,i 1 ) + O(e 2 ). (13) 

Substituting Eqs. lfT2")) and (fT3"|) in Eq. fTTj) and collecting terms at the order e = and 1, we have, 

a^occ + = 0, (14) 

(15) 



where the subscripts refer to derivatives with respect to £ and r\. The solution to Eq. (|T4[) can be taken as 

x a (tv) = A{rj) cos^ + B{? 1 ) sin£, (16) 

where the coefficients A and B are functions of the slow time-scale. Substituting Eq. lfTI))) into Eq. lfTS")) and removing 
the secular terms [l6j |. we have the following coupled differential equations for the slow time-scale, 

A' = -kB +\nA-\nA(A 2 + B 2 ) + \\B, (17) 
B' = kA+-nB -^B(A 2 + B 2 ) + ^XA. (18) 
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Figure 1: Bifurcation diagram of Eqs. (|17[) and (|18[l in the fc-A plane. The line A = 4|fc| along which a Hopf bifurcation occurs 
at the origin. Below this line, in the shaded region, there exists a limit cycle. The dashed line denote a saddle-node bifurcation 
at the origin. 
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Figure 2: Entrainment by the parametric dust-charge fluctuation for a = 0.15, (3 = 0.1, A = 1.0, u>d = 1.0, e = 0.5, where the 
amplitude R of the oscillation is plotted against the period p of the parametric forcing term. The two vertical dashed lines 
at p = 2.8 and 3.6 indicate the region of entrainment (marked as EN). The other behavior is quasi-periodic (marked by QP). 
Only very little hysteresis is observed in the quasi-periodic regions. 

We note that the hyperbolic fixed points of the slow flow correspond to the periodic motion of the original equation 
Eq. (fiTj) i.e. an entrainment and limit cycles of the slow flow correspond to quasi-periodic motions of Eq.(fTT]) 
From simple resonance dynamics it is evident that the entrainment region of Eq. by the parametric forcing term 
should increase as the parametric forcing amplitude A increases, allowing an wider range in the detuning parameter k 
during which the entrainment is observed. Therefore it is worthwhile to study the flow of the slow variables through 
Eqs. fTT]) and (fTg|) as we vary the parameters k and A. 

In Figfl] the bifurcation diagram of the slow flow in the k-X plane is shown. Along the line A = 4\k\, a Hopf 
bifurcation occurs at the origin, below which, in the shaded region, a limit cycle appears. As A falls below 4|fc|, the 
fixed point at the origin becomes an unstable spiral [complex eigenvalues of the linearized Jacobian of Eqs. (fTT|) and 
(fl8|)] from an unstable node. Along the line A = 4\/l + 4/c 2 , saddle node bifurcation occurs for sufficiently high A when 
the origin becomes a saddle point from an unstable node after two other saddle points coalesce at the origin. From 
what we have observed, one can conclude that entrainment by parametric forcing occurs above the shaded region of 
Figfl] In the shaded region, we expect a quasi-periodic behavior. These results are confirmed from the numerical 
solutions of Eq. (fTT]) , which are shown in FigfS] The entrainment region for A = 1 can be calculated from Figfl] as 
—0.25 < k < +0.25, which in terms of period of Eq. pT]) for the parameters of Figfl] is given by 2.79 < p < 3.59 which 
very well agrees with the numerical results (see FigfSJ. 

As we increase e, the expansion parameter, the prediction from Figfl] however agrees less and less as the system 
makes a transition to chaos. 
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IV. STABILITY AND BIFURCATIONS 



The stability of a periodic orbit can be effectively studied using Floquet theory [19|, [20j . In this section, we briefly 
review the Floquet theory with reference to Eq. (|10p . We first write the second order dynamical equation Eq. (|I Op as 
two first order equations, 

x = y, (19) 
y = (a — (3x 2 )y — w d x{\ — e\ cosvi). (20) 

We have already shown in the previous section that a frequency-locked (entrained) phase with a single periodic limit 
cycle can exist for Eqs. (|19I20|) . In this section, we are going to study the stability of this periodic orbit and bifurcations 
leading to the onset of chaos. 

The Poincare map of an initial point zq — (xo,yo) on the periodic orbit (limit cycle) can be obtained by sampling 
the orbit points z n at discrete time interval t = t n , n = 1,2,3, .. .. So the transformation which successively maps 
the Poincare section P(z) is z n+ \ = P(z n ). The linear stability of a g-periodic orbit with P q (zo) = Zq can now be 
determined from the linearized map given by the matrix DP q of P q at an orbit point zo, where P q is the g-times 
iterated Poincare map. The linearized matrix M — DP q can be obtained by integrating the linearized equations 
corresponding to Eqs. (|19l20p for small perturbations along the g-periodic orbit [l9l ]. 

Assume that z*(t) = z*(t+ q) is a point lying on the g-periodic limit cycle of Eqs. (|19|20p . We perturb the orbit 
with a small perturbation 5z = (Sx, Sy) and linearize Eqs. (|19l20p about the closed orbit, 

^ A*) (t), m = (? i) , (2D 



5 v J w V<W \fx f y j ix , y) = {x *, y *) 

where J(t) is the q-periodic linearized Jacobian. The partial derivatives f x , y are given by, 

f x = -2f3xy-cj 2 d {l-e\cosut), f v = (a-(3x 2 ). (22) 

We now assume that W(t) — [wi (t) , W2 (t)] is a fundamental solution matrix with W(0) = I [lj| . The general solution 
of the g-periodic system, Eq. pTjl , is then given by, 

(SS)^KS)' 

We then substitute Eq. ([23l) into Eq. ([2T|) to obtain the initial value problem, 

W(t) = J(t)W(t), W(0) = I, (24) 

where W(q) is the linearized map DP q (zo). So, the matrix DP q can, in principle, be obtained from numerical 
integration of Eq. (|24p over period q. However, the numerical procedure is not very straight forward and requires 
sophisticated techniques to determine the exact form of the matrix DP q which is very sensitive to initial conditions 
(x*,y*). The matrix DP q is known as the monodromy matrix for the g-periodic orbit and the eigenvalues of this 
monodromy matrix, popularly known as the Floquet multipliers [19, 20], indicate the stability of the g-periodic orbit. 
Therefore, the values of the Floquet multipliers have to be determined with considerable precision for understanding 
the true nature of the stability of the nonlinear system. The characteristic equation of the linearized map M = DP q 
is given by, 

C 2 - < + A = 0, (25) 
where the eigenvalues £i 2 are the Floquet multipliers and r = tr(M), A = det(M). The determinant A is given by 

EMU 

A = e So f(J)«« = e (a-/3«* a )9 ) (2 6 ) 

(27) 



- f tr( J) dt (mod 
1 Jo XI 



We know from Floquet theory that the periodic orbit is stable only if the pair of Floquet multipliers lie inside the 
unit circle. The bifurcations of the periodic orbit occur on the unit circle. From the expression for A, Eq. (|26p . it can 
be seen that we do have bifurcations depending on a balancing of the parameters a, /3, and the periodic orbit, which is 
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Figure 3: Stability diagram of Eqs. (|19l20fl in the a-e plane. The period doubling bifurcations occur along the lines denoted by 
the points 'x'. The points denoted by a 'o' along the dashed line indicate pitchfork bifurcations. The chaotic regime lies above 
the two lines. As we can see that some pitchfork bifurcations are followed by a period doubling cascade. The other parameters 
are /3 — 0.1 and v — ujd — A = 1.0. 



Table I: Scaling of period doubling cascades. 
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determined by the parameter e, the magnitude of the parametric driving force. In all probability, the unstable region 
should lie in the region of large a. We numerically determine the Floquet multipliers for a range of periodic orbit in 
the parameter space (a, e) with /3 = 0.1 and ujd = A = 1.0 and the resultant stability diagram is shown in Figj3] In 
FigEl period doubling bifurcations occur along the lines denoted by a ' x ' sign and subcritical pitchfork bifurcations 
occur along the dashed line denoted by a 'o'. The line joining the 'x' points in the figure denotes the accumulation 
points or the limiting points of the period doubling bifurcations before transition to chaos. The chaotic regime lies 
above these lines. As the two lines seem to intersect, when extended, we see that some of the period doubling cascades 
are preceded by pitchfork bifurcations. 

As the pair of Floquet multipliers decreases through —1 at the period doubling points, the g-periodic orbit loses its 
stability to jump to a 2g-periodic stable limit cycle. So, in the observed parameter regime of Fig[3l all period doubling- 
bifurcations are supercritical [l6| . In case of the pitchfork bifurcations, the Floquet multipliers increases through 
+1 and are subcritical as there arc no stable limit cycles after the bifurcations and the system becomes aperiodic 
[l6|. In FigsGUa) and (b), two successive period doubling orbits are shown. In FigsJ3fa) and (b), we have shown the 
bifurcation diagrams with the bifurcation parameters as a and e (for details, please see the captions in the figure). 
These bifurcation diagrams are obtained with the help of AUTO [22j as a part of the XPPAUT package [23| . 



A. Scaling of the period doubling cascades 



It is interesting to investigate the scaling behaviour of the period doubling cascades in light of the scaling of 
the period doubling sequences in 1-D maps. As usually observed in any period doubling cascades, the bifurcation 
parameter, which in our case are a and e, converges geometrically to a limiting value with the convergence ratio 
approaching a unique value, analogous to Feigenbaum number in case of 1-D maps [H, |24|- In Table]]] we have 
listed the values of the bifurcation parameters as these converge to their respective limiting values. The ratio of this 
convergence 5k, expressed as, 

6k= A A k ~ Ak -\ lhn4 = <5, (28) 




Figure 5: Bifurcation diagrams for the period doubling sequences. In (a), the bifurcation parameter is a with e = 1.5 and in 
(b), the bifurcation parameter is e with a = 0.15. All other parameters are same as in Fig[3] The open boxes in the diagrams 
are unstable orbits and the solid lines indicate stable periodic orbits (limit cycles) with period doubling bifurcations. Blow-up 
regions of the second period double cascade is shown in the insets. 



approaches a unique limiting value S as the bifurcation parameter Ak — > const., with k denoting each successive 
period doubling point. The value of 8 characterizes the scaling property of the bifurcation, which agrees well with the 
Feigenbaum number 4.669 . . . for 1-D maps || [LH. |24|]. 



V. TRANSITION TO CHAOS 



In this section, we study the chaotic behavior of Eq.fTU]). The single most prominent feature of chaos is its sensitive 
dependence on initial conditions, which is measured by the Lyapunov characteristic exponents (LCEs) (25j. These 
exponents are invariant global indicators of the non-linear system. For a continuous dynamical system described by 
a set of autonomous ordinary differential equations, the number of LCEs is equal to the dimension of the system. By 
discretizing the temporal dimension as At, the LCEs can be defined as, 

^=^}^oAi Hs * {MN)] > (29) 

where 

N 

M N = ]Je J ^ At , (30) 

4=0 

In the above relations, Si are the singular values of the matrix Mjy and J is the Jacobian. Numerically, the number 
N denotes the number of integration steps of length At. Here, we employ a numerical algorithm, based on the Wolf's 
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Figure 6: (a) Plot of the maximal Lyapunov exponent, cri, for Eq.JTUJ) in the parameter space of a-e. The Lyapunov exponent 
is calculated in the entire parameter space indicated in the figure. The blackened portions indicate a positive and the blank 
portions of the figure indicate a negative Lyapunov exponent. As can be seen from the figure, the fractal nature of the chaos 
is indicated by repetitive appearance of the whole figure in the smaller regions close to the horizontal axis. It also seems that 
the entire figure is only a part of a large figure covering the entire a-e plane, (b) The orbit diagram for Eq. (|31fl . where the 
successive local maxima are plotted against the bifurcation parameter, for a — 0.15. Other parameters are same as in (a). 



well known method to calculate the LCEs 26]. Note that, Wolf's original method is not the best approach to calculate 
the LCEs and it can be modified to contain the non- uniformity-factors of the LCEs [27| • An important necessary step 
in Wolf's algorithm is to re-orthogonalize of the set of vectors, which is carried out, usually, through Gram-Schmid 
orthogonalization procedure. In our numerical routine too, we have used the Gram-Schmid orthogonalization to 
calculate the LCEs. In particular, we have calculated the LCEs for the following autonomous system, 

x = y, 

V = (a- (3x 2 )y - u>%x(l- e\u), 

u = u(l - u 2 - v 2 ) - 2nv/T, 1 ' 

v = v(l - u 2 - v 2 ) + 2-ku/T, 

where v = 2t:/T. Note that, we have made the original non-autonomous equation, Eq. (|10[) . a system of autonomous 
equations by introducing the last two equations of Eqs. pTj) . which have stable and unique solutions, 

u(t) = cos(2vrf/T), v(t) = sin(27rf/T), (32) 

for the initial values [m(0), f(0)] = [1,0]. 

In FiglHJa), we have plotted the maximal LCE for the Eqs. (f3"Tj) . in the parameter space of a-e for the range shown. 
The parameters are f3 = 0.1, u 2 , = 1, and A = 1. The period of the orbit is chosen as T — 6.5 which corresponds 
to v = 0.96664. The blackened portions of the figure indicate a positive Lyapunov exponent indicating chaos. In 
all other places the maximal Lyapunov exponent is negative, signifying stable oscillations of the system. We can see 
the fractal behavior of the chaos from the figure. It also seems that the whole figure is only a part of a larger figure 
covering the entire domain of the a-e plane. In Fig[71 the maximal LCE, <n, is plotted along with a blow up of the 
region of period doubling cascades. The period doubling points are marked with 'arrows' in Fig[7£b). 

We have constructed an orbit diagram [l6[ [Fig|Bfb)] for the system, Eq. (j3"Tj) , where we have plotted the successive 
local maxima for the variable x of the oscillation against the bifurcation parameter e. The period doubling bifurcations 
occur near 1.9 and 2.2 before the system becomes chaotic. Note the appearance of 5-period window in the region 
e = 2.5 and 3.0 after which the system becomes chaotic again. The chaotic orbit along with the time evolution of the 
system in the chaotic regime and the Poincare map are shown in Fig [5] 



VI. CONCLUSION 



In this work, we have carried out a detailed investigation of the stability, bifurcation leading to chaos for the vdPM 
system arising out dust-charge fluctuation in a a dusty plasma. We have shown that the system can be highly chaotic 
depending the chosen parameters and not as restrictive as has been pointed out by Saitou and Honzawa [lj. In fact, 
an wide range of chaotic region exists for the parametric driving strength (eA) as low as 0.05, as shown by the plot of 
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Figure 7: (a) Plot of the maximal Lyapunov exponent <ti against the bifurcation parameter e for a = 0.15. Other parameters 
are same is in Fig[6] (b) A blow up of the previous figure is shown where the period doubling points are marked with arrows. 
Note that after each period doubling bifurcation, the system transits to a stable limit cycle (characterized by large negative a\ 
before getting unstable for the next period doubling point when o~\ becomes close to zero, (c) The Maximal LCE in a region 
of pitchfork bifurcation. The first arrow (up) indicates the pitchfork bifurcation after which the system transits to a chaotic 
regime before becoming periodic (denoted by the negative LCE) and becomes chaotic again as a increases, following a period 
doubling cascade, denoted by the second arrow (down). 

the maximal Lyapunov exponent o\. It has also been found that the system can be completely deterministic in the 
middle of two chaotic regions and exhibit quasi-periodicity, as shown by the orbit diagram where a 5-period window 
appears in the middle of a chaotic region. We have further shown that when the parametric forcing term related to 
dust-charge fluctuation is small, away from the chaotic regime, the system can be driven in a frequency-locked state 
when a harmonic resonance of 2:1 takes place between the driving frequency and the fundamental frequency of the 
system. We have found that in most of the cases, the system transits to chaos through a cascade of period doubling 
bifurcations and the scaling of the period doubling cascades closely agrees to that of 1-D maps [24j . 
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